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1. INTRODUCTION 

Decomposition is the basis for modeling complex seasonal time series. A time series is said to have 
a complex seasonal pattern when it has trend and multiple seasonal patterns, perhaps with non-integer period. 
Recently, this kind of series has been intriguing researchers to study and develop methodologies to improve 
the forecast accuracy. Soares and Medeiros [1] have discussed TLSAR (two level seasonal autoregressive) 
model that consists of two parts, deterministic and stochastic component. This method is a simple statistical 
model that combines linear trend, trigonometry and the autoregressive model. Meanwhile, [2-7] developed 
an innovation state space exponential smoothing model. Taylor [2] added the second seasonal component 
to the Holt-Winter exponential smoothing, namely HWT, and in the following year, [6] presented BATS 
(Box-Cox transform, ARMA errors, trend, and seasonal component) and TBATS (trigonometric BATS) to 
reduce the parameters of HWT model. Different from other previous methods, TBATS has accommodated 
the fractional seasonality in the data. Still, it unable to take into account the time varying component. 
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Meanwhile, other researchers [8-11] considered SSA method in complex time series forecasting. 
The development of the SSA method was inspired by [12-15]. Much of subsequent works refer to their 
works. Later, [16-23] continued to develop SSA as a tool of time series forecasting. Vautard et al. [16] 
reviewed the capability of SSA in extracting information from short and noisy time series without any prior 
knowledge about the data. Yiou et al. [17] regarded SSA as a method to identifying and extracting oscillatory 
component from the original series. Vahabie et al. [18] combined SSA with autoregressive (AR) model to 
predict the electricity demand. In the following year, [19] and [20] reviewed and demonstrated the success of 
SSA in economic and financial time series forecasting. Golyandina and Korobeynikov [21] described 
the RSSA package and showed how the methodology of SSA can be implemented for analyzing, estimating 
parameter and forecasting. Hassani et al. [22] studied the effect of outliers and showed that the vector SSA 
forecasts were more robust than the recurrent SSA forecast. Khan and Poskitt [23] discussed the properties of 
forecast obtained by SSA in theory and practice. Recently, [24] and [25] implemented SSA method to 
decompose a series into several components and modeled the components by the Box-Jenkins method. 
However, some data exhibit a deterministic behavior rather than stochastic. Thus, there will be parts that 
have not been accommodated by the Box-Jenkins model. Recently, [26] to combine SSA-LRF with NN 
(neural network) and fuzzy time series to forecast hourly electricity load in Indonesia. The success of NN in 
dealing the nonlinearity can also be seen in [27-31]. Moreover, the method proposed by [24-26] did not 
take into account the seasonal with non-integer period, time varying amplitude, and calendar variation. 
Hence, some modifications are needed to accommodate the features of the series instead. 

This work proposes two SSA-based hybrid model, namely TLSNN (two level seasonal neural 
networks) and TLCSNN (two level complex seasonal neural networks) model, modifications of TLSAR 
model. TLSAR consists of deterministic and stochastic component. The deterministic is expressed as 
a combination of linear, harmonic, and dummy variables to take into account the linear trend, seasonal, and 
holiday influences. Meanwhile, the two proposed models provide flexibility in determining deterministic 
component, especially in describing nonlinear trends. Moreover, one of the proposed models can capture 
the oscillation pattern with time varying amplitude. In this study, the two models are constructed by 
considering the results of SSA decomposition. We take advantage of the capability of SSA in decomposing 
time series so that we get several separable and interpretable components including the information on how 
much each component contributes to the series. Further, SSA method makes the identifying process of 
the component patterns easier. Our methodology is as simple as that in [1]. The main difference is that we 
take into account the time-varying components of the series that are not included in TLSAR model or other 
methods mentioned in the literature. Further, we can show that our methodology can improve the multi-step 
forecasting accuracy. 


2. THE PROPOSED METHOD 
2.1. Mathematical definition for the proposed SSA-based hybrid model 

The two proposed models are called the two-level seasonal neural networks (TLSNN) and the two-level 
complex seasonal neural networks (TLCSNN). Both models are modifications of the TLSAR model. 
The higher order of polynomial trend and the time-varying sinusoidal function are included in 
the deterministic component to capture a more complex pattern in the time series. In this case, the results 
from the SSA decomposition make the identification and determination of the proper deterministic function 
for each component of TLSNN and TLCSNN easier. Generally, TLSNN and TLCSNN are expressed in 


Y; = 54 + Zt (1) 
Where Y, is observation at time f, S,; is the deterministic component and Z; is the stochastic component. 
The difference between the two is in the form of S+. Meanwhile, Z+, is approximated by neural network (NN). 
These two models are explained further in Definition | and Definition 2. 


Definition 1. A time series {Y,,¢ = 1, ..., N} with the sample size N and which has a complex pattern follows 
TLSNN if the deterministic component in (1) can be presented as 


c= ye aitt + X, [aon cos(2ttht /T) + Bon sin(2mht/T)] + X}? 8xDgt (2) 


i=0 


and the stochastic component is 


Z=) b Oy CZ) + er (3) 
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where 


2 


E Ga) 1 (4) 


Ninp 7 
1—exp(-2 dj CjZt-j) 


The notation @;, Aon, Pon» Ôk for i = 1, ..., npon A = 1, ... Np and k = 1,2, ...,Np are the unknown parameters 
of the model, T is the seasonal period, n, is the number of harmonics, and np is the number of dummy 
variables included in the model. D,; notates the effect of the kth category at time t. 

Definition 2. A time series {Y;,¢ = 1,...,N} with a complex pattern follows TLCSNN if the deterministic 
component in (1) can be presented as 


Npol . 
= yo Gey yy ao eats cos(w,t) + Bygt? sin(w,t)| 


5 
+ Drilon cos(2rht/T) + bon sin(2rht/T)| + ye Ôk Dkt ~ 


and the stochastic component is as in (3) and (4). 

The notation Qj, Qyg,Byq,%on Pon Ox for i= 1,...,Nnpob Y = 1, nosih = 1, .. Nh q = 1, ..., npo, and 
k = 1,2, ...,Nnp are the unknown parameters of the model, T is the seasonal period, np is the number of 
harmonics, and np is the number of dummy variables included in the model. D;; notates the effect of the kth 
category at time t. 

The TLSNN and TLCSNN models accommodate, not only trend linear, but also quadratic or 
another polynomial with higher order to capture other trend behaviors. The differences between TLSNN and 
TLCSNN are in their harmonic terms. We add the time-varying amplitude sinusoid component, do estimate 
its parameters separately using iterative ordinary least square to estimate the parameter of the amplitude 
modulated sinusoid model (see [32, 33]) and then combine its best function with the deterministic component 
of TLSAR model. Meanwhile, the parameters of the deterministic model in TLSNN are estimated 
simultaneously using householder transformation. This estimation method was discussed detail in [34-36]. 


2.2. Modeling Strategy for the two proposed SSA -based hybrid model 
2.2.1. TLSNN model 
The steps of TLSNN modeling procedure are explained below. 
Level 0: the series is decomposed by SSA method into several separable components 
Level 1: obtaining the deterministic model (2) based on the results of SSA decomposition 
a. identifying the trend function and the dummy variables to accommodate calendar effects. 
b. estimating the parameters of the deterministic model by using least square based on 
householder method [34-36] 
c. defining the fitted value based on the deterministic model, $, 
Level 2: obtaining the stochastic model 
a. determining the stochastic component, Z, = Y, — S; 
b. testing the linearity in the series by using Terasvirta method [37] 
c. modeling by NN when the linearity test detecting the nonlinearity relationship in the data 
evaluating the accuracy of the forecast values 


2.2.2 TLCSNN model 

The procedure of TLCSNN modeling are explained in the following steps. 

Level 0: the series is decomposed by SSA method into several separable components 

Level 1: obtaining the deterministic model (5) based on the resuts of SSA decomposition 
a. identifying and estimating each component function 

— estimating parameters of the trend function by ordinary least square (OLS) method 

— estimating parameters of the oscillatory function by iterative OLS [33] 

— defining the appropriate function for representing the irregular component. This function 
may involve Fourier function with the dummy variables. The number of harmonics can be 
obtained based on AIC value. 

b. defining the fitted value based on the deterministic model, $, 
Level 2: obtaining stochastic model 
a. determining the stochastic component, Z; = Y, — $, 
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b. testing the linearity in the series by using Terasvirta method [37] 
c. modeling by NN when the linearity test detecting the nonlinearity relationship in the data 
evaluating the accuracy of the forecast values 


3. RESEARCH METHOD 
The SSA-based hybrid model is defined by the method that combines deterministic and stochastic 
component, where the deterministic component is obtained based on the SSA decomposition result. 


3.1. The benchmark model 

The TLSAR is used as the benchmark model in this research. TLSAR model is presented as 
a combination of a trend linear and trigonometric function of time and the autoregressive model. The TLSAR 
model is expressed as the combination of deterministic and stochastic component as in (1), 
where the deterministic component is written as 


Se = ao + at + E alaik cos(2mtht /P) + Bon sin(2mht/P)] + X;? 8kDrt- (6) 


The notation ao, Q4, Qon: Pon, Ox for h = 1,...,n, and k = 1, 2,...,Np are the unknown parameters of 
the model, P is the seasonal period, np, is the number of harmonics, and np is the number of dummy 
variables included in the model. The number of harmonics nj, in (6) is determined by minimizing the Akaike 
Information Criteria (AIC), 


AIC = nln(6*) + 2M, (7) 


go Gai f l l l 
Where ô=- £r (Y; — Y,)% is the estimate of the residual variance. M is the number of parameters 


fitted in the model, and n is the number of residuals that can be calculated from the series (see [34, 38, 39]). 
The dummy variable D;, identifies the effect of the kth category at time t, for example, the kth hours of 
the day to Y; or the kth days of the weeks to Y,. This value of D;; 1s one if the kth category is observed at 
time t and zero for otherwise. In the original TLSAR model, Z; 1s the stochastic component that follows the 
AR model while in this study, Z; follows the SAR model. The autoregressive model of order p, denoted by 
AR(p) is given by 


Pp (BDZ: = er, (8) 


Where ¢,(B) = 1 — ppB — + — @,B? and e; is zero mean white noise process. This model can be 
applied to the series in which the preceding values and the random shock influence its present value. 
The process needs to be stationary so that the roots of ¢,(B)=0 must lie outside the unit circle. 
When the seasonal part is stochastic and cannot be handled by the deterministic model, (8) can be 
extended as 


Py (B)Pp(B*)Z, = êt, (9) 


Where p(B) = 1 — pB — = — pB”, pp(B°) = 1 — p1(B°) — +: — pp (B®), s is the seasonal 
period and P is the order of the seasonal component. The model presented in (9) denoted by 
SAR(p, 0,0)(P, 0,0), 


3.2. SSA method in decomposing a time series 

SSA method decomposes an original time series into trend, several oscillatory components, and a 
noise. The major references related to this method can be found in [11, 40-43]. This method consists of four 
steps, embedding, singular value decomposition (SVD), grouping and diagonal averaging. The time series 
{Y,,t = 1,2,...,N} with the sample size N is transformed into a Hankel matrix Y with L- dimension vectors, 
that is 


Y Y baa Ve 
y= "2 "3 o "+1 (10) 
Yo Year ot Yy 
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Where L is a window length and K = N — L + 1. This transformation process is called by 
the embedding step. The embedding step is then followed by calculating YY’ and applying SVD to determine 
its eigentriples Wii, U;,,V;) for i = 1, 2, ..., L. The A; represent the i-th eigenvalue of matrix YY’ and JA 
denotes the i-th singular value of matrix Y. The U; and V; are the left and the right singular values of matrix 
Y. SVD of the Hankel matrix in (10) can be written as 


Y=Y¥,+ Y+-°+Y,y (11) 


where d is the rank of matrix Y and Y; = Jij U,V; for j = 1, 2, ..., d. The contribution of the matrix Y; in 
(11) can be determined by calculating the ratio A;/ yi Aj. The matrices (Y1, Y2, ..., Yq) are arranged in m 
groups, where m < d such that (11) can be written as 


Y = Y, + Yi + eee + Yii (12) 


In this case, each matrix Y;, in (12) may be reconstructed from one or more eigentriples (4/4;, Uj, Vj) and the 
contribution of this matrix can be measured by the ratio )ijé,, A;/ pan Àj. 


Each matrix Y;, is then retransformed into a time series by the diagonal averaging procedure. 
The diagonal averaging is a Hankelization process that transforms the elements of matrix Y;,, say y;j into the 
t-th observation of the i-th (i = 1, 2, ..., m) component series by averaging y;; for all i, j such thati + j = 
t + 2. Thus the original series {Y,,t = 1, 2,..., N} can be represented as the sum of m series, 


(= (13) 


where y® is the series reconstructed from the matrix Y;,. 


The separability of all component series, i i = 1,..,m, in (13) can be measured by calculating w- 
correlation values, which are in the range of [-1,1] and visualize them in the form of a w-correlation matrix 
diagram with different level of color from the brightest (white means no correlation or separable) to 
the darkest (black means perfect correlation). The strongest separability is represented by the white color and 
the weakest separability is shown by black color (see [40]). Other alternatives to achieve separability to 
obtain a more parsimonious forecasting model were proposed by [44] and [45]. 


3.3. Diagnostic checking and evaluating the performance of forecast accuracy 

Model diagnostic checking is needed to check whether the model assumption that is zero mean 
white noise residual series is satisfied. To see whether the residuals are uncorrelated random shock with zero 
mean and constant variance, we can examine both the plot of residuals and its sample ACF. We can also 
verify the results by using the Ljung-Box test [46, 47]. In this work, we use RMSE and MAPE to measure 
the performance of the forecast accuracy. 


RMSE = E n (Y, — Î,)2, (14) 


1 


Yt—-Y, 
MAPE = -Xiz | 


Yt 











(15) 


where n is the number of observations involved in the calculation. Y, and Y, are the real value at time t and 
the forecast value at time t, respectively. Y is the mean of the in sample data. Equation (14) and (15) are used 
for evaluating the training data and the testing data. 


4. RESULTS AND DISCUSSION 

The well-known monthly accidental deaths in the USA and the daily load series of Jawa-Bali for 
hour 1.00 am to hour 4.00 am were used in this empirical study. In this section, we show how the two 
proposed method implemented to the complex series and compare the results with TLSAR as the benchmark. 


4.1. Monthly accidental deaths in USA 


The monthly accidental deaths series depicted in Figure 1(a) is quite popular and has been discussed 
in several studies. Brockwell and Davis [48] have used this data as an example of the application of 
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the ARIMA, Holt-Winter seasonal, ARAR model, while [11] discussed the application of the SSA-LRF 
model. The data are divided into two parts, the training and the testing data set. We used the data from 
January 1973 to December 1978 as the training data and January 1979 to June 1979 as the testing data. 
As in [11], window length ZL used in this study is 24, proportional to the period. Based on the w-matrix 
correlation, the time series can be decomposed into 4 components (see, Figure 1(b)), namely trend, two oscillation 


components, and noise. 
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Figure 1. The monthly accidental deaths series from January 1973 to June 1979, 
(a) The original series, (b) Its decomposition 


In TLCSNN model, the quadratic function is the most suitable for modeling trend components, 
while the first and the second oscillation components can be approximated by quadratic amplitude modulated 
sinusoidal function and linear amplitude modulated sinusoidal function, respectively. Since Terasvirta test 
indicates that there is a nonlinearity relationship in the stochastic component, NN must be more appropriate 
to model the stochastic component than the AR or other linear models. The best NN is determined based on 
the RMSE value. The number of input variables are chosen between 6 or 12 regarding the seasonal periods 
while the number of hidden units are chosen between | to 10 so that the network produces the smallest 
RMSE value with random residuals. Based on the experimental results, NN with 6 input units and 8 hidden 
units, denoted NN (6-8-1) is the best fit NN model for the stochastic component. 

Further, in TLSNN model, we also choose the quadratic trend in the deterministic component and 
NN(6-8-1) for the stochastic component. In this case, the number of harmonic is 5 with period 11.98. 
The comparison results of RMSE and MAPE for the three models are summarized in Table |. It can be seen 
that TLSNN and TLCSNN produce smaller RMSE and MAPE than TLSAR both for the training and testing 
data. TLSNN and TLCSNN respectively reduce the RMSE of TLSAR by 85.26% and 82.44% for 
the training data and by 82.73% and 95.00% for the testing data. Meanwhile, sequentially TLSNN and 


TLCSNN reduce the MAPE of TLSAR by 86.69% and 87.15% for the training data and by 82.52% and 95.58% 
for the testing data. 


Table 1. Comparison of RMSE and MAPE obtained by TLSAR, TLSNN, and TLCSNN 
model for the death series 


Model Training Testing 
RMSE MAPE RMSE MAPE 
TLSAR (h=2, P =12, SARIMA(0,1,1)(0,1,1)12 270.4575 2.4997 334.7116 3.1314 
TLSNN, Myo = 2, Nn =5, NN(6-8-1) 39.8595 0.3327 57.7906 0.5474 
TLCSNN, npoi = 2, Nos = 2, nn =0, NN(6-8-1) 47.4913 0.3213 16.7392 0.1384 
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4.2. Daily electricity load of Jawa-Bali in the specific hours 

The electricity load of Jawa-Bali, in megawatt per hour (MWh), obtained from a power utility 
company in Indonesia, Jawa-Bali control center (JCC), that were observed hourly are considered in this 
work. However, the data discussed is only the electricity load at 01.00 am, 02.00 am, 03.00 am, and 04.00 
am. These data have a unique pattern compared to other hours due to the influence of the habits of Indonesian 
citizens in the month of Ramadan. These four periods are considered to represent complex seasonal patterns, 
which have seasonal patterns, trends, and calendar influences. The data pattern of time series 01.00, 02.00, 
03.00, and 04.00 for the period of January 1, 2009, to December 31, 2011, is presented in Figure 2. 
Each series consists of 1095 observations. It is divided into two parts, the first 1088 observations are the 
estimation sample (the training data) and the last 7 observations (a week) are the test sample. 
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Figure 2. Daily load series from 1 January 2009 to 31 December 2011 for, (a) 01.00, (b) 02.00, 
(c) 03.00, and (d) 04.00 


There was a noticeable surge in several days and repeated from year to year. Higher surges occur at 
03.00 and 04.00. This surge occurred during the sahur hours in the month of Ramadan and eventually 
dropped dramatically before and after the Eid al-Fitr. Each series exhibits trend and seasonal. The demand for 
electricity that increases from time to time is related to population growth and society needs. This causes 
a positive trend in the series and can be modeled by the function of time. Due to the similarity results in 
the component patterns of SSA decomposition for 01.00 to 04.00, we only report 01.00 in detail. Nevertheles, 
the overall results will be dicussed. The results of SSA decomposition for time series 04.00 is shown in 
Figures 3(a) and 3(b). 

Based on Figure 3(c), the dummy variables in the TLSAR, TLSNN, and TLCSNN models can be 
described in Table 2. In this case, the electricity load tends to fall on the next day (the day after a holiday). 
Especially for the month of Ramadan, the electricity load increases. The k dummy variable {k = 1, 2,..., 10} 
is one if time t shows the event as described in Table 2 and has a value of O for the other. All the models 
described in Table 3 have harmonic components with the number of harmonics n, = 3 and dummy variables, 
as stated in Table 2. The results of the comparison of RMSE and MAPE are presented in Table 3. 
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Figure 3. The components of series 04.00 obtained by SSA decomposition with L = 490, 


(a) Trend component, (b) Oscillation component, (c) Irregular component 


Table 2. Types of days presented as dummy variables in the TLSAR, TLSNN, and TLCSNN model 


Description 

The second day of Ramadan up to five days before the first Eid Al-Fitr 

The two days of Eid Al-Fitr and the two days after the second Eid Al-Fitr 

The third day after the second Eid Al-Fitr 

A day before the first Eid Al-Fitr = and the fourth day after the second Eid Al-Fitr 

The second day before the first Eid Al-Fitr, and the fifth and sixth day after the second Eid Al-Fitr 

The third day before the first Eid Al-Fitr, a day before the first Eid Al-Fitr, and the seventh day after the second 

Eid Al-Fitr 

A day after the independent day, new year, and holidays related new year. 

Eid al Adha and a day after Eid Al Adha 

A day after Chrismas, a day after Waisya, a day after Nyepi, a day after Imlek, a day after the good Friday, a day after the 

ascension day of Jesus, a day after New Hijr. = 

A day after presidential and legislative elections 

Table 3. RMSE and MAPE calculated from the three models for the training and testing data 

Method RMSE MAPE 
Training Testing Training Testing ( 25-31 December 2011) 
1 2 3 4 5 6 7 

TLSAR 258.1585 604.5183 1.5269 3.0608 2.2952 2.6805 2.7826 2.7215 3.3189 3.8516 
TLSNN 225.0217 210.5599 1.3292 1.7816 1.1340 1.6498 1.4100 1.4120 1.3194 1.3108 
TLCSNN 219.6763 223.3640 1.2854 1.9569 1.8204 2.0931 1.6147 1.3931 1.2624 1.3573 
TLSAR 258.4460 598.5510 1.5406 1.9241 1.5728 2.2643 2.5783 2.7039 3.3912 3.8528 
TLSNN 226.9396 291.9397 1.3586 1.5159 1.8073 1.8611 1.9852 1.9391 1.6678 1.9340 
TLCSNN 226.5588 284.2629 1.3412 0.7976 0.8453 1.3011 1.9835 2.0730 1.7419 1.5428 
TLSAR 258.0119 594.6462 1.5426 2.6900 1.8375 2.4459 3.0090 3.0498 3.6538 3.9846 
TLSNN 228.1649 364.7621 1.3492 3.0535 1.6601 2.0184 2.9127 2.3784 2.2030 1.9297 
TLCSNN 216.0240 282.1139 1.2868 1.1977 1.1545 2.0739 2.3221 2.1597 1.8226 1.5671 
TLSAR 166.6101 674.3093 1.0052 3.2022 1.6647 2.1011 2.8500 3.1066 3.7152 4.3406 
TLSNN 230.2454 431.8189 1.3081 2.6571 1.6434 2.0005 2.7185 2.4458 2.7858 2.7745 
TLCSNN — 218.3380 229.3032 1.2888 1.6704 1.6563 1.3914 1.9303 1.7635 1.4701 — 1.2627 


Based on Table 3, it can be concluded that the TLSNN and TLCSNN models provide a more 
accurate forecasting value than the TLSAR model. The modulated oscillation component in the TLCSNN 
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model increases the forecast accuracy of the TLSNN model, especially at 03.00 and 04.00. The TLSNN 
model is able to reduce the value of the RMSE of TLSAR between 35.96% and 65.17% for the testing data, 
while the TLCSNN model is able to reduce the RMSE of TLSAR between 52.51% and 65.99% for 
the testing data. For the case of electricity load forecasting at 01.00, 02.00, 03.00, and at 04.00, the TLSNN 
and TLCSNN models are more acceptable than the TLSAR model. Furthermore, the results of the TLSNN 
model forecasting for 02.00 and the TLCSNN model for 04.00 provide MAPE values that meet PLN 
standards, which are below 2%. 


5. CONCLUSION 

In this study, we proposed two SSA-based hybrid forecasting methods for complex seasonal time 
series and provided five data series to show the application examples. The proposed methods named TLSNN 
and TLCSNN were developed from TLSAR by considering components of SSA decomposition results. 
The proposed hybrid models consist of deterministic and stochastic components. The deterministic 
components of TLSNN and TLCSNN models constructed based on the SSA decomposition results 
accommodated not only trend linear but also quadratic or another polynomial with higher order to capture 
other trend behaviors. The differences between TLSNN and TLCSNN were in their harmonic term. 
Parameters of the deterministic model in TLSNN were estimated simultaneously using householder 
transformation. Meanwhile in TLCSNN, we added the time-varying amplitude sinusoid components to 
the trend and harmonic or stationary sinusoids. The parameters for each functions were then estimated 
separately. Finally, the stochastic component of TLSNN and TLCSNN were modeled by NN to overcome 
nonlinearity relationship problem in the data. Being modified, the models tend to produce better performance 
for the multi-steps forecast values than TLSAR for the case of the monthly accidental death series and Jawa- 
Bali load series. The findings of this study have several implications for future practice. Another possible 
future research would be to develop the methods for multivariate problems. 
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